Theories, models, simulations: a computational challenge 



G.C. ROSSP 

Dipartimento di Fisica, Universita di Roma Tor Vergata, 
INFN, Sezione di Roma Tor Vergata 
Via delta Ricerca Scientifica, 00133 Roma, Italy 

In this talk I would like to illustrate with examples taken from Quantum Field Theory and Biophysics 
how an intelligent exploitation of the unprecedented power of today's computers could led not only 
to the solution of pivotal problems in the theory of Strong Interactions, but also to the emergence of 
new lines of interdisciplinary research, while at the same time pushing the limits of modeling to the 
realm of living systems. 



Prologue 

The somewhat schematic partition of the last century natural science into separated fields 
of research, which were essentially identified with mathematics, physics and biology, is 
nowadays becoming less and less rigid, leading to large areas of overlapping interests. 

The fundamental reason for the former de facto separation was the enormous amount 
of accumulated knowledge in each of the three areas, which resulted in an increasing, and 
at the end unsurmountable, degree of specialization for people working at the front-end 
of their research field. 

Two facts have been drastically changing the situation. One was the growing evidence 
that methods and ideas developed in one research area could be fruitfully exported to 
other, even distant, fields of investigation. The second, is a related one and has to do with 
the sharp increase of the available computing resources (in terms of CPU-time, memory 
and storing capacity), which is making algorithms and general computational strategies 
immediately ready for use to researchers working in different areas. 

In my opinion this last fact is of particular relevance in today's spectacular progress 
of science, because it has allowed to imagine and attack problems that were considered 
impossibly difficult only a few years ago. New, flexible and adaptive computational tools 
that can be of general help to many scientific disciplines are being implemented under the 
pressure of the challenges posed on the one hand by the developments of pure science 
and technology and on the other by the fast expanding needs of our modern societies 
(think to weather forecasting, stock market "surveillance", power plant control systems, 
distributed information network management, etc.). 

Taking an example of this trend from a field which is nearer to the scientific inter- 
ests of our community, it is interesting to remark that one of the most extraordinary and 
somewhat unexpected outcome of the long lasting interplay between Statistical Mechan- 
ics and the theory of Strong Interactions in its lattice formulation (lattice QCD - LQCD) 
was the decision taken within the community of theoretical physicists to build "dedicated 
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machines" with parallel architecture . The aim of these machines was to provide a 
tool capable of efficiently dealing with the extremely hard computational task of extract- 
ing useful physical information from the simulation of QCD, when the latter is seen as 
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a statistical system of interacting "coloured spins" living on the sites of a (Euclidean 
space-time) lattice with gauge fields sitting on the links. 

Numerical and conceptual tools developed in Statistical Mechanics and in Theoretical 

Chemistrjl^J^ immediately found applications in LQCD, and vice-versa ideas and nu mer- 

ical techniques invented in LQCD were fed back in simulations of statistical systems0^ 
as well as in the study of the more complicated situations that appear when systems of 



biological interest are modeledl^El 

The second half of the 80's was marked by a breakthrough in the theory of disor- 
dered systems that turned out to have a significant impact in numerous emerging fields 
of investigation. The replica approach was extended to spin glass systems and the notion 
of replica symmetry breaking was proposed as an explanation for the occurrence of the 

glassy phase transition ^ . In this context a new and more precise notion of complexity 
has emerged, suggested by the phenomenology of spin glass, that rather soon appeared 
to be of great relevance in the apparently distant problem of constructing mathematically 
sensible models of biosystems. 

In fact, there is an intriguing analogy between the mathematical structur e of spin 

glasses'^ and certain approaches to the problem of modeling protein folding E PI^^ I 
Here the relation between the two fields is in physical terms less direct than in the case 
of Quantum Field Theory and Statistical Mechanics mentioned above and most impor- 
tantly in the case of both spin glasses and proteins mathematical computationability is 
intrinsically limited by the complexity of the models one is considering. Despite these 
difficulties, a lot has been learned about protein structure from approaches inspired by 
the theory of disordered systems and, vice-versa, ideas taken from biology have spurred 
new strategies aimed at dealing with hard computational problems (NP-complete prob- 

from a novel point of view. 
Indeed, the very recent discovery that the "typical" (not the worst) NP-complete prob- 

lerrJ^ (examples of NP-complete problems are the i^T-SAT problems^^) may be (almost 
always) solved with polynomial algorithms (like the cavity method, or the "survey in- 
spired decimation") seems to suggest and make us hope that similar methods could be 
developed and used to attack the most challenging among the theoretical problems that 
arise in modeling biological macro-molecule interactions, among which we might men- 
tion protein folding and aggregation, protein-protein and protein-DNA recognition, etc. 

I cannot end this brief overview without recalling that the most spectacular and suc- 
cessful results of the generalized use of computers in pure and applied research are prob- 
ably to be found in the realm of life sciences. Sequencing the human genome would 

have been impossible without the support of the most advanced computers of the time!^. 
Today the big task is annotation. It is now clear that to gain a really useful understanding 
about the complexity of living systems, we need to record, cross-link and organize in an 
appropriate way the exponentially fast growing amount of biological information that is 
being gathered in experiments. The task is made particularly difficult by the impressive 
variety of data we need to store and correlate. Just to give you some examples of such 
an enormous variety let me recall that understanding biosystems at large will require 
dealing with data that go from the structure of the metabolic networks of biochemical 
reactions taking place in the cell to the description of the series of events by which the 
immunological system responds to an antigen, from the epidemiological and statistical 
information necessary to monitor the progression and the spreading of a disease in a pop- 
ulation to the biochemical characterization of the complicated protein-DNA interactions 



which regulate gene expression and so on. 

The very same development of the micro-array technique, that so much biological in- 
formation is continuously providing, was only possible thanks to the wide-spread avail- 
ability of computers capable to deal with the huge outflow of data of combinatorial chem- 
istry in an efficient, reliable and retrievable way. 

1 Introduction 

Personally I was introduced to the fascinating field of computers and simulations by 
Adriano in 1980, when we were both visiting CERN. It was the exciting time when the 
first attempts to extract physics from numerical simulations of QCD were just starting 
to produce useful results and APE was a new extraordinary scientific and technological 
enterprise. 

Since then the increase of the computational power at disposal to research and every- 
day life has proceeded at a pace that only the most blunt extrapolation of the Moore law^ 
over more than thirty years was able to predict. This exponential explosion has radically 
changed not only the life style of billions of people, but also the way we scientists think 
about science and research. Completely new problems have appeared to be within our 
reach, that only few years ago would have seemed just impossible to attack or even to 
dream of. If appropriately used, computers represent more than a simple tool which 
can increase our ability to answer questions: their enormous potentiality, associated to 
flexibility and adaptability, has opened the way to new adventures that are only limited 
by our fantasy and courage. 

In this talk I would like to try to underline the irreplaceable role of what might be 
called "intelligent computing" in certain domains of physics and biophysics, by illus- 
trating in three significative examples of application, chosen according to my personal 
inclination and competences, how new ideas could be effectively implemented and made 
to work thanks to the power of the available computational means. Two examples are 
taken from the field of Monte Carlo simulations of LQCD. The first has to do with the 
analysis of the gluon sector of QCD (sect. The second with possible ways of solving 
or easing the problems posed by the explicit breaking of chiral symmetry which accord- 
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ing to the Nielsen-Ninomiya theorem ° affects any (ulti-a-)local lattice regulaiization of 
QCD (sect.|51l. In the third example I wish to report on a somewhat innovative approach 
to the study of polymer structure with the methods of Statistical Mechanics (sect.@}. 

2 Gluon operators 

I want to start by discussing two selected topics related to the gluonic sector of QCD 
where "intelligent computing" has been decisive to give support to our understanding of 
certain properties of the Theory of Strong Interactions. I will illustrate the calculation and 

Moore's original statement was the observation made in 196^3 that the number of transistors per square 
inch on integrated circuits (we would more precisely say today the number of transistors that minimizes the 
cost per transistor in a chip) had doubled every year since the integrated circuit was invented. Moore, co- 
founder of Intel, predicted that this trend would continue for the foreseeable future. In subsequent years, 
the pace slowed down a bit, but transistor density has doubled approximately every 18 months, and this is 
the current definition of Moore's Law, which Moore himself has blessed. Most experts, including Moore 
himself, expect Moore's Law to hold for at least another two decades. 



the physical relevance of two quantities: the plaquette expectation value and the topolog- 
ical susceptibility. The first quantity is related to the so-called gluon condensate The 

second is supposed to be responsible^^ for the non- vanishing of the 77' mass in the chiral 
limit (the limit where quark masses are sent to zero). 



2. 1 The plaquette expectation value 

The expectation value of the plaquette, (P) , is an obviously relevant quantity in the study 
of the thermodynamic properties of lattice gauge theories. Besides, it was thought that 

one could extract the F^-gluon condensate of ref. from lattice data if one could subtract 

from the lattice data on (P) its perturbative tail^^ . In this context it was an open question 

to decide whether signs of renormalon effects0l and of what dimension were visible in 
the plaquette perturbative expansion. 

At the time where we (I mean Adriano and me) started to ask ourselves such questions 
there was little experience about perturbative and non-perturbative definition of lattice 
composite operators and even less about the relation between lattice and continuum ex- 
pectation values. Lacking any better strategy, we attacked by brute force the problem of 
defining the F^-operator starting from its definition in terms of the plaquette expectation 
value. We computed the first three terms {i.e. tree-level, order (1-loop) and order g'^ (2- 
loops)) in the perturbative expansion of the plaquette by hand. At that time ours was the 
most difficult perturbative lattice calculation ever attempted. It took us about six months 
of intense work and cross-checking until we could agree on the analytic expression of the 

function that we then had to integrate numerically!^. The result was 
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where Nc is the number of colours and the error in parenthesis comes from the uncertainty 
inherent in the numerical integration. It is amazing (should I say disappointing looking 
back at our effort?) to observe that the clever stochastic methods which are available 

todayl^ allow to compute the perturbative expansion of {P) up to order {g^Y^, with the 

aid of a good 32-node cluster in a few hours . From this recent knowledge indications 
are that a dimension four operator can indeed be seen below the computed perturbative 

tail, if accurate data, like those of ref.122] are used in the analysis. 

To tell the truth we could do something slightly better: by comparing our results to 

the brand new Nc = 2 simulation data just produced in those days by Mike Creutz , we 
could extract the numerical value of the coefficient of the term g^ obtaining an estimate 
which, within its relatively large error, appears to be quite accurate, when compared to 
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the successive explicit perturbative calculations of ref. . 

This story is paradigmatic of the inextricable interplay between technological devel- 
opments and scientific intelligence. Thanks to his scientific creativity Mike Creutz was 
able to exploit at their best the computational possibilities of the time, producing data 
that led us to ask questions whose answer had in turn to wait still a few years before one 
could arrive at the technical improvements and theoretical advances necessary to get a 
full comprehension of the underlying problems. 



2.2 Topological charge density and susceptibility 

Topology is a key concept in gauge theories. According to our understanding of the 
solution of the so-called U(1)a problem a non- vanishing topological susceptibility is re- 
sponsible for providing a mass to the t]' pseudo-scalar meson in the limit where up, down 
and strange quark masses are set to zero. As a result the rj' is not the ninth Goldstone 
boson of chiral symmetry. 

In the limit Nj/Nc — > ^ the mass of the (lightest) flavour singlet pseudo-scalar 



meson is given by the well-known Witten-Veneziano (WV) formula!221 

where Fj^ is the pion decay constant (normalized so that F^^ ~ 94 MeV for A'^^ = 3) and 
A is the "topological susceptibility". A is formally defined by the equation 
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with Q{x) the topological charge density, which in the formal continuum theory has the 
expression 

Qix) = E F^FpAx) ■ (4) 



a=l 



The notation (. . .)|ym in eq- © means that the QQ-correlation function is to be com- 
puted in the pure Yang-Mills theory, i.e. in the absence of quarks. 

The idea that the non-perturbative value of A could be measured from pure gauge 

lattice simulations dates back to the works of ref.l^, where the first attempts to extract 
such a number from numerical data were made. The resulting quantity, though non- 
vanishing and endowed with the correct scaling behaviour, was yielding a value of the t]' 
mass significantly smaller than phenomenologically required. 

The discrepancy was due to the fact that the renormalization effects necessary to 
match lattice and continuum definitions of topological charge density had been com- 
pletely overlooked. This mistake was corrected in the seminal paper of ref.l^, where 
the required renormalization constant was computed to one-loop in perturbation theory. 
Remarkably when the perturbatively normalized and vacuum subtracted simulation data 
for the topological susceptibility were inserted in eq. Q, the agreement between the the- 
oretical calculation and the experimental value of the rj' mass turned out to be rather 
good. 

In my opinion getting an agreement between theory and experiments in this comer 
of the theory is of especially great conceptual importance, because the r]' mass issue is 
one of the few instances where the non-perturbative structure of QCD as a theory for 
Strong Interactions is at stake and can be subjected to a stringent test. For this very good 
reason the Pisa group (led by Adriano) has striven for some time to arrive at an accurate 
and fully non-perturbative definition of the topological objects relevant to this problem. 
Indeed they have been finally able to get a reliable non-perturbative determination of 



'^Nf is the number of light (massless) quark flavours. 



the renormalization constant and subtraction term necessary to construct from simulation 
data the proper definition of A. Tliis was achieved relying on the clever method of cooling 

the gauge configurations to freeze out their perturbative fluctuations)^. 
2.3 Topology in chiral regularizations ofQCD 

The situation of LQCD simulations has radically changed recently owing to the appear- 
ance on the market of exactly chiral fermions '^ ^ i | 32,33 , 34i , jj^g subsequent observa- 
tion that the index theorem holds true as a lattice identityl^^ if fermions obeying the 
Ginsparg-Wilson (GW) relation are employed. In this framework the WV formula 

can be given a rigorous non-perturbative status . In fact, after identifying the unrenor- 
malized operator which represents the topological charge density on the lattice as the 
one suggested by the flavour singlet Ward-Takahashi identities of the GW-regularized 
theory, one can prove that eqs. (|3j and ^ are valid on the lattice with no need for any 
renormalization or subtraction. 

The trouble with this approach is that simulations where the definition of A suggested 
by GW fermions is employed are fairly expensive, although rather nice results have been 

recently obtained for it . Adriano's recent idea in this context is surprisingly simple 
and effective: it consists in making use of the GW-inspired definition of topological 
charge density only to the extent the latter is needed to determine the non-perturbative 



normalization constant of the more standard gluon definitionE21 i.e. only to measure the 
topological charge of a configuration. The interest of this strategy is obvious: it allows 
to get an accurately normalized topological charge density without having to pay a much 
too high computational price. 

3 Waiting for a fully chiral simulation of LQCD 

The next generation of computers may allow LQCD simulations with exactly chirally 
invariant fermions, i.e. fermions obeying the GW-condition . In the meanti me a vi- 
able alternative'^ could be to employ maximally twisted Wilson fermions 4 2|43 | 44 | 
possibly accompanied with a judicious choice of the pure gauge action. Preliminary 

quenched as well as unquenched numerical results in this direction are quite en- 
couraging. They confirm the theoretical expectation that correlators are 0(a) improved 
and that simulations require computational times t hat ar e of the same order of magnitude 
as for plain Wilson fermions (see, however, sect. 13.21 for some word of caution). Ex- 
trapolation of the present trends makes us confident that the overall computational power 
allocated in Europe to maximally twisted lattice QCD (Mtm-LQCD) simulations can 
match the CPU-time needed for a study of the full theory in physically realistic condi- 
tions, i.e. on a (3 fm)^ x6 fm lattice with a pion mass of about 250 MeV. The computation 
requires an estimated power of the order of 10 Teraflop*year'^. Optimistically one may 
hope to get the first useful results in a little more than one year from now. In view of 
this remarkable and fortunate situation I think it might be worth reviewing the theoretical 

structure and the properties of Mtm-LQCD as developed in refs.l ^^ l ^^ l ^^ l ^ l 
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3.1 A cheap proposal 



Soon after noticing that to avoid exceptional configurations in Wilson fermion simula- 
tions one should introduce quarks in flavour pairs and have the Wilson term rotated with 
respect to the quark mass term by an axial rotation in iso-spin space, it was realized that 
an especially useful choice for that angle is to set it at its maximal value, \uj\ = 7r/2, 
because in such a situation 0(a) (actually 0(a^'^+^), /c > 0) improvement of physical 
quantities is automatic with no need to introduce the "clover term"'^in the action. 



It was then shown irp2lthat the nice improvement properties enjoyed by Mtm-LQCD, 

which were derived for pairs of mass degenerate quarks in^^ can be immediately ex- 
tended to the more interesting case of non-degenerate quarks, without loosing the posi- 
tivity of the corresponding fermion determinant. The last property is obviously crucial 
if one wants to be able to set up workable Monte Carlo-like simulation algorithms for 
LQCD. 

With the above ingredients and exploiting the flexibility offered by the freedom of 
regularizing different valence flavours with different values of the Wilson parameter, it 

was shown inl^ that it is possible to construct a hybrid theory, where sea quarks are 
introduced as pairs of non-degenerate particles and valence quarks are regularized as 

Osterwalder-Seiler fermions, such that no "wrong chirality" mixing affects the 
computation of the matrix elements of the CP-conserving AS* = 1,2 effective weak 
Hamiltonian. Of course the same result would hold if GW fermions were used as valence 
quarks. Absence of wrong chirality mixing makes Mtm-LQCD a more appealing regu- 
larization of QCD than the one offered by the use of standard (clover) Wilson fermions. 

From what we said above about improvement, it turns out that Mtm-LQCD correlators 
that are not trivially vanishing in the continuum limit can be affected by lattice artifacts 

described by a Symanzik expansionESl with only even powers of a. Among these terms 
there are lattice contributions which tend to become large as the quark mass is lowered. 
They originate from the breaking of parity and iso-spin induced by the presence of the 
twisted Wilson term in the action. These effects have been discussed both in chiral per- 
turbation theoryl^ES as well as in the language of the Symanzik expansion*^ where 
they appear as terms of the form (a/mg)^^, k > 1. The general conclusion of the theo- 
retical analysis is that such lattice artifacts can be reduced to a numerically tolerable level 

(precisely down to order a^(a^/mg)'^"^, /c > 1) if the clover temJ^is introduced in the 
actioJ^ (with its non-perturbatively determined csw c oefficiently^ or, alternatively, if 
the critical mass is chosen in some "optimal way"I^E2l44. Actually it turns out'^^'that, 
at least up to 0(a) included, the optimal critical mass coincides with the critical mass one 
would get from the vanishing of the pion mass (or the PCAC mass) within the standard 
Wilson fermion regularization. 

The previous discussion about chirally enhanced discretization artifacts affecting Mtm- 
LQCD correlators is rather important because it shows that the strong (order of magni- 
tude) inequality 

> «^QCD > (5) 

invoked in ref.l^ in order to have the phase of the chiral vacuum driven by the quark 
mass term and not by the (twisted) Wilson term, can be relaxed to the more favourable 



relation 



(6) 



before large cutoff effects are possibly met when the quark mass is lowered at fixed a. 
The bound Q is fairly weak as it permits simulations in a region of quark masses that 
correspond to rather light pions (with masses around 200 MeV for typical present-day 
lattice spacings). 



3.2 Where is the catch? 

All this sounds good, perhaps too good to be true. So the natural question to ask is: is 
there a catch in the twisted mass approach to LQCD and where is it? 

To tell the truth there is one little catch. It has to do with the observatioJ^^^^^that 
at too coarse lattice spacing meta-stabilities are seen to affect unquenched data -J which 
prevent their extrapolation to the chiral limit. Such meta-stabilities are the consequence 
of the explicit breaking of chiral symmetry induced by the presence of the Wilson term 
in the action. They appear at sufficiently low quark mass when the latter is progressively 
lowered at fixed a and cause the statistical system one is dealing with not to reach equi- 
librium. For recent reviews on these and related questions and an updated assessment of 

the present status of quenched and unquenched Mtm-LQCD simulations see ref.l^. 

A safe way-out of these difficulties is obviously to work at sufficiently small lattice 
spacing: something which, however, may turn out to be computationally too expensive. 
Actually there is another, more clever solution to the existence of meta-stable phases 

(and to the other dangerous flavour breaking effects 1521 that in this regime plague the 
theory) which will work even on coarse lattices. It consists in tuning the pure gauge 
action so as to set to zero (in the chiral limit) the matrix element of the dimension six 
operator of the Symanzik low energy action of LQCD taken between pion states with 
vanishing three-momentum, i.e. the quantity C2 oc (7r(0)|£6|^(0))- It can be shown that 
this particular matrix elem ents controls the magnitude of all the unwanted cutoff effects 

described abovel^Q^' ^^l^^i making them to vanish as soon as C2 = 0. 

This strategy has been already partially implemented by working with a gluon action 

other that the standard plaquette action '^^ One finds that meta-stabilities are avoided 
in this way as soon as a < 0.1 fm and for pion masses down to (at least) 300 MeV. 



4 Structural properties of polymer chains 

Lacking at the moment in most cases mesoscopic, functionally useful, descriptions of 
biological systems, theoretical models aimed at understanding the dynamic and/or the 
thermodynamic properties of molecular aggregates of biological interest are based on 
a detailed atomistic description of the compound. The physico-chemical behaviour of 
the resulting model and its compatibility with the available experimental information is 
then investigated by numerical simulations. The deterministic approach of Molecular 
Dynamics (MD) or the stochastic methods of Monte Carlo type^^^J^llSJ are employed 

either classically or with quantum corrections injected d la Car-Parrinello . 

The need for a numerical approach appears to be even stronger if the problem of 
predicting the folded configuration of a protein, solely from the knowledge of its lin- 
ear amino-acidic composition, is considered. The interest of investigating the folding 



problem rests on the experimental observation that the biological functionality of a pro- 
tein crucially depends on the nature of its folded configuration'^. Misfolding is, in 
fact, known to lead to malfunctioning and in certain cases to severe pathologies, such as 

Creutzfeld-Jacobs disease and human variant of BSE 64, Alzheimer disease^, cystic 

fibrosis and probably also to other neuro-degenerative processes. 

Understanding the nature of folding is expected to be a formidable task: already the 
classical problem of finding the absolute minimum of the free energy of atomistic models 
of long polymer chains has a computational complexity which bears close resemblance 

to that of instances belonging to the class of problems technically called NP-complete^! 
Furthermore the problem may not have a unique solution: the recent studies on misfold- 
ing induced deseases have shown that proteins may live in more than one (meta-)stable 

state. It is remarkable that the simple model of ref . ^ ^ can yield some understanding for 
this behaviour. 

In the following sections I shall discuss merits and limitations of some interesting re- 
search lines and computational strategies that have been recently put forward to deal with 
the problem of folding or, more modestly, with the problem of predicting the structure of 
a polymeric chain from its chemical composition. For a review of approaches of different 

nature see, for instance, ref.^^ 



4.1 State of the art 

The study of Statistical Mechanics of polymers, i.e. long chains consisting of monomers 
of specific nature, is becoming more and more important in chemical technologies and 

biological applications. Polymers like proteins'^, nucleic acids'^, polysaccharides'^^ 

and synthetic materials'^ display features that strongly depend on their detailed physico- 
chemical properties like, for instance, the degree of flexibility of certain chemical bonds, 
the charge density at monomer atoms, the structure of the hydrogen bond network be- 
tween monomers either close or far away in the sequence, and so on. 

As we said, the enormous complications associated, even within classical physics, 
with the atomistic description of the specific interaction among the elementary compo- 
nents of the polymer can only be handled by numerical simulations. Clever algorithms 
have been devised to explore the configuration space available to the system and differ- 
ent types of ensembles have been invented and numerically implemented, starting from 

molecular dynamics (MD) and Monte Carlo (MC) methods ^1. As is well known, MD 
and MC simulations explore the micro-canonical ensemble and the canonical ensemble 
of the system, respectively. Other kinds of ensembles, which may be collectively indi- 
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cated by the name of generalized ensembles , have also been introduced and employed 
for the study of thermodynamic properties at equilibrium. In principle, under standard 
ergodicity assumptions, all these ensembles should yield equivalent physical informa- 
tion. Indeed the use of generalized ensembles within MD and MC simulation strategies 
resulted in a rather powerful approach capable of predicting the statistical properties of 
atomistic models of fluids and other compounds of chemical and/or biological interest in 
a wide range of temperatures and order parameter valuesP^I. 

The crucial limitation that is encountered in numerical simulations of systems with 
many relevant degrees of freedom is related to the inadequate and strongly biased sam- 
pling of the configurational space occurring when the temperature is lower than the criti- 



cal temperature of the model. By "critical temperature" we generically mean the temper- 
ature above which the system is in the disordered phase. Below the critical temperature 
the system remains trapped in local minima, within energy barriers that are rarely (or 
never) overtaken by thermal fluctuations. 

Many strategies have been proposed over the years ai med at trying to overcome this 
difficulty. Among them we may recall simulated annealin^^^Ell stochastic tunnelin^25J 



and many variants of genetic algorithms'^. The problem with these approaches is that 
they do not always permit the calculation of statistical averages in well defined ensem- 
bles {i.e. in the ensembles that are statistically representative of the desired experimental 
conditions). 

Vice- versa, algorithms designed to access ensembles of the multi-canonical type, i.e. 
the kind of generalized ensembles that exploit the information on the (potential) energy 

density of states^^lSl are well ass essed and made rather effective if used in conjunction 

with the replica-exchange method^ ^ l ^ ' ^^ i In many interesting instances it is possible to 
computationally monitor order parameters of geometrically constrained molecular mod- 
els of polymers in a fairly large temperature range. 

However, even within these approaches problems arise when all the many degrees of 
freedom of reahstic models, including the high-frequency vibration modes, are taken into 

account ^^J, as it is necessary to do in order to treat condensed phases and explicit sol- 
vents. In fact, large variations of the potential energy are observed associated with such 
stiff terms in the Hamiltonian even for tiny configurational changes. Such large poten- 
tial energy changes cause very low acceptance in the exchange of temperatures between 
replicas and lack of convergence in the multi-canonical weight computationl^ES 



4.2 A proposal for a new approach 

The main lesson one learns from the previous discussion is that energy is not the best 
variable to label configurations because on the one hand configurations that are only 
slightly different in their atomic spatial arrangement may have largely different potential 
energies and on the other configurations with similar energy can be structurally very 
different. This is the main reason why, within the standard multi-canonical approach, 
introducing the temperature through the modulation of the energy density of states by the 
Boltzmann factor does not yield sufficiently satisfactory results, as soon as the number 
of degrees of freedom of the polymer is too large and/or the temperature is above the 
order-disorder phase transition. 

In order to overcome this type of problems it was proposed in'^ to work in a gen- 
eralized ensemble where configurations are generated according to the density of states 
associated to some configurational quantity (or some set of configurational quantities), 
rather than energy. These configurational quantities can be the mean value of some bond 
or dihedral angle along the polymer chain, the value of the a-helicity of the polymer, the 
head-to-tail distance of the chain or any other variable which may serve to characterize 
the geometrical structure of the system. 

A problem with this approach may be considered the fact that it is not clear how 
to introduce the notion of temperature, because the fundamental statistico-mechanical 
relation between the energy of the system and its temperature is put at stake. Actually, 
for the purpose of studying, say, biopolymers, which after all are neither isolated systems, 
nor do they work at equilibrium, this state of affairs is not really a problem and can be 



dealt with along the lines described below in sect. 14. 31 (see also the Appendix). 

In brief the idea of ref.'^ is to start by working in the micro-canonical ensemble 
associated to some configurational variable, A, rather than energy, and then pass to the 
associated canonical ensemble by the "constrained maximal entropy method" (CMEM) 
imposing that A assumes some preassigned value, a. The latter can be either taken from 
experiments or can be known from some exact theoretical calculations in particularly 
simple models. The scheme can be extended in an obvious way to the case of more than 
one configurational variable. 

The passage from the standard m icro -canonical/canonical ensembles to the configu- 
rational ensembles introduced in ref.'^is schematically illustrated in the steps 1. to 4. 
outlined below (in the formulae that follow we generically indicate with r the whole set 
of variables necessary to describe the degrees of freedom of the system). 

1. Make the replacement 

U{r) A{r) 

oj^{E) = J dr5{E - U{r)) ujA{a) = J dr5{a - A{r)) 

2. From seeds at random temperatures (see Appendix) collect configurations accord- 
ing to a Metropolis test with multi-canonical weight 

obtaining the configurational distribution (which may or may not be further elabo- 
rated) 

Pu{r) ^ pA{r) 

3. Determine the best configurational distribution, P, satisfying the constraint 

{U)=JdrU{r)Pu{r-J)=E ^ {A) = J drA{r)PA{r;~X) = d , (7) 
using the CMEM, which yields 

^c/(r;/3) = ^Pc7(r)e-'5f^M ^ A) = ^P^We-^^W 

ZuiP) = J drPu{r)e-^^''^^ Za{X) = J drPA{r)e-~^^^'''^ 

with the Lagrange multiplier implicitly fixed by Q 
P = P{E) A = A(a). 

4. The expectation value of F = F{r) is computed by means of the re-weighting 
formula 



where 

E, = U{ri), F, = F{ri) ^ ai = A{r,) , Fi = Fin) , 
and Nconf is the number of collected configurations. 



(8) 



Roughly speaking we may say that the computational strategy displayed in the left col- 
umn is fine for energy related quantities, but not so much for structural quantities. On the 
contrary the new strategy outlined in the right column is expected to work appreciably 
well for structural quantities, but not as well for energy related quantities. To appropri- 
ately deal with them temperature must be brought back on stage. 

4.3 Introducing temperature 

Introducing the notion of temperature for a complex (fully flexible) system, like a poly- 
mer, is a delicate issue, because of the observation we already made that configurations 
only slightly different in their atomic spatial arrangement may have largely different (po- 
tential) energies. Consequently, as it turns out, it becomes more and more difficult to get 
the correct (Boltzmannian) energy distribution of the total available energy among the 
many degrees of freedom of the system as the temperature increases (despite the fact that 
at high temperature overcoming energy barriers may become easier). 

Actually, in the scheme we have just discussed there is room for the introduction 
of a sensible notion of temperature. This is done in two separate, but complementary 
steps. Temperature can be injected in the configurational probability distribution. Pa, 
if we know how the expectation values of the configurational variable we have chosen 
to fix {i.e. a) depend on T. This dependence will in turn induce a T dependence in the 
values of the Lagrange multipliers that are obtained by solving the constraint equation Q. 
Through eq. (ISjl this dependence is then passed to the expectation value of any other 
configurational quantity one wishes to compute. It is important to remark that in a similar 
way dependence upon other environmental parameters can be introduced in the study of 
the physico-chemical properties of the system. 

The temperature dependence induced through the method described above is not 
enough, however, to produce the correct T behaviour of quantities that require an ac- 
curate thermalization of all the degrees of freedom of the system for their calculation. 
Examples of such quantities are the moments of the (potential) energy distribution. In 
these cases a local, extra thermalization step has to be carried out. This can be accom- 
plished in the following way. Starting from each one of the recorded configurations, 
one performs a number of hybrid MC steps with velocities extracted from a Maxwell- 
Boltzmann distribution at the desired temperature. At the end of each MD block of moves 
configurations are subjected to a standard Metropolis test with acceptance/rejection prob- 
ability given by exp{—[5H), where H is the total (kinetic plus potential) energy of the 
system. In this way configurations are smoothly thermalized at the desired temperature 
and can be used to compute the ensemble averages (|SJl. 

4.4 An application to oligopeptides 

A first step in the study of protein folding properties can be the determination of the 
local propensity of the amino-acid chain to form a-helix, /3-sheet or other more or less 
structured arrangements. As an application of the considerations previously illustrated in 

this section, I would like to briefly report on the interesting example considered in ref.^ll 
where the propensity to form a-helix structures of two simple oligopeptides, viz. Glyi2 
(a chain formed by 12 Glycine amino-acids) and Alai2 (a chain formed by 12 Alanine 
amino-acids), was studied. The result of comparing data from the simulations of the two 
oligopeptides was that the order (folded) to disorder (unfolded) critical temperature is 
lower for Glyi2 than for Alai2, implying that the propensity of Glyi2 to form a-helix 
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Figure 1: Difference of average potential energy between disordered (Ua = 0) and ordered (na = 12) states 
for Alai2 (squares) and Glyi2 (circles) as function of temperature. 



is lower than for Alai2. The conclusion Tc{G\yi2) < Tc(Alai2), which is in agreement 
with experimental evidence, follows from Fig. 1 by identifying Tc as the temperature 
at which the energy difference between the ordered and the disordered phase equals the 
equipartion energy. For completeness I give in the Appendix some detail on how the data 
of Fig. 1 are obtained following the strategy described in sect. 14.21 

I wish to end this section with two remarks. First of all the result one gets is rather 
robust and does not depend on the precise definition of critical temperature one decides to 
employ, as Alai2 data all lie higher than Glyi2 data. Secondly and, most importantly, the 
possibility of determining such a non-trivial physico-chemical property of these systems 
should be regarded as a methodologically rather remarkable result because it is based on 
a first principle computation. It only relies, in fact, on the topology of the systems and 
the detailed properties of the atomistic model we took to describe the force field of the 
two oligopeptides. 

5 Conclusions 

It was my intention in this talk to convey to you the idea that the enormous computational 
power we have today at our disposal has been a crucial ingredient for the spectacular ad- 
vances we have been witnessing in many research areas as well as in every-day-life appli- 
cations and technological developments. I tried to do so by showing in a few significative 
examples, belonging to my field of competence, how the cultural layout upon which new 
ideas and methods have emerged have been appreciably influenced by the easy access to 
large-scale computing facilities and vice-versa. 

In recent years a sort of mini-cultural revolution has been taking place. The radically 
reductionistic paradigm, that so successful had proved to be in our quest for the fun- 
damental laws of micro-physics, is appearing not to be fully adequate to deal with the 
challenge posed by the physics of, say, dynamical (non-linear) systems or the conceptual 



problems of modeling biological objects. Chances are that in these emerging fields of 
investigation notions like chaos or complexity will be going to play a central role. These 
notions have rapidly evolved from the initial physico-mathematical frameworks where 
they have been first introduced (non-linear dynamics and the physics of disordered sys- 
tems). They have grown to the status of conceptual interpretative schemes under the 
stimulating pressure of the many successes they have led to in difficult numerical prob- 
lems and the beneficial effect of the vast diversification of their field of application. 

Outlook 

Let me conclude with a personal note. I started my career as a student of Bruno Touschek 
who, as you certainly know, was the inventor of e+ — e~ colliders and the leader of the 
group of scientists that in Frascati built the first working storage ring, christened AdA, 
for "Anello di Accumulazione". A large fraction of theoreticians and students in Rome 
and Frascati were at that time (from 1966 to 1970) busy with computing cross-sections 
for all sorts of e+ — e~ processes. Computers were not based on transistors or chips, 
but on electronics tubes, and people were still busy with punching cards. So you were 
expected to cany on your calculations as much analytically as possible and only at the 
very end come up with some sensible approximation that one could work out numerically. 
Nevertheless Touschek was fond of electronic computers. He had clear in his mind their 
enormous potential in many strategic applications and he immediately suggested that 
computers should be used to simulate statistical systems with the idea of checking theory 
against actual numerical data. 

To a large extent this the same competent and inspired vision I found in Adriano's 
attitude towards research. One can identify a clear and consistent line of development 
in Adriano's scientific activity which, starting from his innovative studies on the role of 
topology in gluon-dynamics at finite and vanishing temperature, has naturally brought 
him in his more recent papers to attack the most difficult problem of all in QCD, the 
problem of understanding the mechanism underlying colour confinement (for a recent 

review see, for instance, ref-ISHl). 

I'm pretty sure that for many years to come Adriano's enthusiasm for research will 
still be a stimulating example for all of us: his ideas and intellectual ingenuity have had 
an enormous impact in lattice QCD, and more broadly in the whole field of high energy 
physics. 

Collaborating with Adriano was a privilege for me which has strongly influenced my 
approach to physics and shaped my scientific interests in research. I wish to thank him 
for that, but most of all for his invaluable and sincere friendship. 
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Appendix 

The definition of ordered and disordered thermodynamical states is given within the con- 
text of the strategy described in sect. 0]by going through the following steps. 



1) Configurations (seeds) are initially generated by sequential MD moves of fixed 
length by taking as starting system coordinates the coordinates of the last stored configu- 
ration and, as starting particle velocities, vector components extracted from a Maxwell- 
Boltzmann distribution at a random temperature, uniformly chosen at each MD step 
within zero and a high-temperature limit (1000 K in the case at hand). This procedure 
can be proved to obey the detailed balance principle and generates a time independent 
(stationary) conditional probability, Pc- Although unknown, Pc is perfectly well defined 

and gives rise to an acceptable probability distribution, 'P^'^\ 

2) As a configurational variable relevant for the problem one is studying the average 
molecular a-helicity, Na, is naturally taken. Na is defined as the number of amino-acids 
with the two dihedral angles, C{i - 1)-N(i)-Ca(i)-C(i) {(pi) and N{i)-Ca{i)-C{i)-N{i + 1) 
(V'i) within appropriately chosen bounds, which were taken to be 260° < 4>i < 320°, 
i = 2, . . . , 12, and 293° < < 353°, i = 1, . . . , 11. 

3) The initial probability distribution, V'^'^\ is improved by multi-canonical iterations 
leading from V'^^^ to 'p('^+^) by generating configurations that are accepted or rejected 
according to a Metropolis test based on the current a-helicity number of states of the 

system, a;j^^(ncj). The iterative procedure is stopped when some stability criterion is 

fulfilled and the last probability distribution, "Pw^, is recorded. 

4) For each oligopeptide the configurational probability distributions corresponding 
to the ordered and disordered phases are constructed from Vnc, by the CMEM, imposing 
the constraint = 12 or = 0, respectively. 

5) At this point the two resulting probability distributions, VN^i^i^a = 12)) and 
{'^{^a = 0)), are thermalized at a set of temperatures ranging from 50 K to 650 K, 

in steps of 50 K. From the latter the data points of Fig. 1 are obtained. 
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